//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//Application to increase the performance of the clusters molecule characterization method implemented by Giovanni Cincilla (IRB-PCB-UB - Barcelona 2008)  //
//This method try to use the initial approximation of one matrix of data to decrease the time and number of iteration. The objetive is to use the complete //
//resource of computation, and parallel the independent task in the algorithms. This programs use the master-worker paradigm, and mpi functions.           // 
//Author: Yandi Naranjo Basalo - IRB - Barcelona - 2009                                                                                                    //
//Date of Last Change: 08-10-2009                                                                                                                          // 
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <iostream>
#include <fstream>

#include <mpi.h>

#include "SmileStore.h"
#include "SmileProfile.h"
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//Constanst to generate the deep of ancestor
#define NANCESTOR    4
#define NITERATION   25
#define PERCENTCONV  0.95
#define LINGOLENGTH  4

//Constanst to create the cluster in algorimths
#define	ARRAYSIZE    60000
#define CLUST_SIZE   3000
#define ITERATIONS   1
#define SHARED       1
#define BUFFCICLE    10000
#define LENGHTBLINE  40
#define NONZERO      0
#define MEDIAN       -0.946667
#define USE_MEDIAN   0
#define USE_DISTANCE 1
#define USE_PREFEREN 1
#define USE_MATLAB   1

#define MAXITS_1     500
#define MAXITS_2     1000
#define CONVITS_1    50
#define CONVITS_2    100
#define DAMFACT_1    0.5
#define DAMFACT_2    0.9

//Name of APA output file 
#define APSTOPFILE   ap.stop
#define IDXOUTFILE   idx.txt
#define NETSIM2FILE  netsim2.txt
#define NETSIMFILE   netsim.txt
#define PREFFILE     pref.txt
#define SUMMARYFILE  summary.txt
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
using namespace std;
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double obtainLingoSim(char *lingoString1, char *lingoString2) {

  long length1 = strlen(lingoString1), length2 = strlen(lingoString2);
  char tmpLingo1[length1 + 1], tmpLingo2[length2 + 1];
  long i, j, k, z;

  //Change the Cl by L and Br by R in the two lingo sequences
  long higher = ((length1 > length2)?length1:length2);

  for (i = 0, j = 0, k = 0; i < higher; i++) {
    if (i < length1) {
      if      (lingoString1[i] == 'l') tmpLingo1[j - 1] = 'L';
      else if (lingoString1[i] == 'r') tmpLingo1[j - 1] = 'R';
      else                             tmpLingo1[j++]   = lingoString1[i];
    }
    if (i < length2) {
      if      (lingoString2[i] == 'l') tmpLingo2[k - 1] = 'L';
      else if (lingoString2[i] == 'r') tmpLingo2[k - 1] = 'R';
      else                             tmpLingo2[k++]   = lingoString2[i];
    }
  }
  tmpLingo1[j] = '\0'; tmpLingo2[k] = '\0';

  j -= (LENGTH_LINGO - 1); k -= (LENGTH_LINGO - 1);
  long   L = j + k; //Total Lingos in two sequences
  double similarityValue = 0, Nai, Nbi;

  for (i = 0; i < j + k; i++) {

    char   compare[5];
    bool   analize = true;

    Nai = 0, Nbi = 0;

    if (i < j) {

      strncpy(compare, tmpLingo1 + i, 4);
      compare[4] = '\0';

      long zHigher = ((j > k)?j:k);

      for (z = 0; z < zHigher; z++) {
        if (z < j && !strncmp(compare, tmpLingo1 + z, 4)) {
          if (z < i) {
            analize = false;
            break;
          }
          Nai++;
        }
        if (z < k && !strncmp(compare, tmpLingo2 + z, 4)) {
          Nbi++;
        }
      }
    }
    else {

      strncpy(compare, tmpLingo2 + i - j, 4);
      compare[4] = '\0';

      long zHigher = ((j > k)?j:k);

      for (z = 0; z < zHigher; z++) {
        if (z < j && !strncmp(compare, tmpLingo1 + z, 4)) {
          analize = false;
          break;
        }
        if (z < k && !strncmp(compare, tmpLingo2 + z, 4)) {
          if (z < i - j) {
            analize = false;
            break;
          }
          Nbi++;
        }
      }
    }

    if (analize == true) similarityValue += (((Nai - Nbi) >= 0)?(1 - ((Nai - Nbi)/(Nai + Nbi))):(1 - ((Nbi - Nai)/(Nai + Nbi))));
    else                 L--;
  }

  return similarityValue/L;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int verifyAPA(char *similaritiesFolder) {

  fstream  filein;
  int      logAPA = 0;
  char    *stateFile = new char[strlen(similaritiesFolder) + 20];
  
  sprintf(stateFile, "%s/ap.stop\0", similaritiesFolder);
  filein.open(stateFile, ios::in);
  filein >> logAPA;
  filein.close();
  if (stateFile) delete[] stateFile;

  return logAPA;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int verifyAPAMatLab(char *similaritiesFolder) {

  fstream  filein;
  char    *line      = new char[1001];
  bool     condition = true;

  filein.open(similaritiesFolder, fstream::in);
  while (!filein.eof()) {
    filein.getline(line, 1000);
    if (strstr(line, "Warning: Algorithm did not converge")) {
        condition = false;
        break;
    }
  }
  filein.close();
  delete[] line;

  return condition;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int loadIdxAPA(char *similaritiesFolder, IdxCluster *IdxArray) {

  fstream  filein;
  int      position = 0;
  char    *stateFile = new char[strlen(similaritiesFolder) + 20];

  sprintf(stateFile, "%s/idx.txt\0", similaritiesFolder);
  filein.open(stateFile, ios::in);
  if (stateFile) delete[] stateFile;
  
  if (!filein.is_open()) { if (stateFile) delete[] stateFile; return 0; }
  
  char *line =  new char[101];
  
  while (!filein.eof())  { 
    filein.getline(line, 100);
    if (strlen(line)) {
      int idxValue = atoi(line);
      IdxArray[position++].IdxValue(idxValue);
    }
  }

  delete[] line;
  
  filein.close();
  return 1;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int loadIdxAPAMatLab(char *similaritiesFolder, IdxCluster *IdxArray) {

  fstream  filein;
  int      position = 0;
  char    *line     =  new char[100];

  filein.open(similaritiesFolder, fstream::in);
  if (!filein.is_open()) return 0;

  while (!filein.eof())  {
    filein.getline(line, 100);
    if (strlen(line)) {
      int idxValue = atof(line);
      IdxArray[position++].IdxValue(idxValue);
    }
  }

  filein.close();
  delete[] line;
  
  return 1;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int computeAPAMatLab(char *similaritiesFolder, char *preferenceFolder, double median, IdxCluster *IdxArray, int length, int posPart) {

  //Load the path of matlab binary programs
  int     convergence = true;

  char   *matlabPATH  = new char[1001];
  fstream fileMatLabConf;
  
  fileMatLabConf.open("matlab.conf", fstream::in);
  fileMatLabConf.getline(matlabPATH, 1000);
  fileMatLabConf.close();

  //Create the temporal file with specific name for each cluster
  char *fileLSIMclusterName      = new char[100];
  char *fileLSIMclusterResult    = new char[100];
  char *fileLSIMclusterStatistic = new char[100];

  sprintf(fileLSIMclusterName,      "LSIMclusterTrial%03d.csh\0", posPart);
  sprintf(fileLSIMclusterResult,    "LSIMclusterResult%03d\0",    posPart);
  sprintf(fileLSIMclusterStatistic, "LSIMclusterStatistic%03d\0", posPart);

  //Dinamic Construction of the script to execute the apcluster with matlab
  fstream fileSIMCluster;
  fileSIMCluster.open(fileLSIMclusterName, fstream::out);

  fileSIMCluster << "unsetenv DISPLAY" << endl;
  sprintf(fileLSIMclusterName, "LSIMclusterTrial%03d\0", posPart);
  fileSIMCluster << matlabPATH << " -r " << fileLSIMclusterName << " >&! matlab" << posPart << ".log << EOF " << endl;
  fileSIMCluster << "plot(1:10)" << endl;
  fileSIMCluster << "print file" << endl;
  fileSIMCluster << "exit" << endl;
  fileSIMCluster << "EOF" << endl;
  fileSIMCluster.close();

  //Dinamic construction of matlab code
  sprintf(fileLSIMclusterName, "LSIMclusterTrial%03d.m\0", posPart);
  
  fstream fileSIMClusterM;
  fileSIMClusterM.open(fileLSIMclusterName, fstream::out);

  fileSIMClusterM << "load " << similaritiesFolder << endl;
  if (USE_PREFEREN) {
    fileSIMClusterM << "load " << preferenceFolder << endl;
    fileSIMClusterM << "[idx,netsim,dpsim,expref]=apcluster(" << similaritiesFolder << "," << preferenceFolder << ",'maxits'," << MAXITS_1 << ",'convits'," << CONVITS_1 << ",'dampfact'," << DAMFACT_1 << ");" << endl;
  }
  else {
    fileSIMClusterM << "[idx,netsim,dpsim,expref]=apcluster(" << similaritiesFolder << "," << median << ",'maxits'," << MAXITS_1 << ",'convits'," << CONVITS_1 << ",'dampfact'," << DAMFACT_1 << ");" << endl;
  }
  fileSIMClusterM << "save('" << fileLSIMclusterResult << "','idx','-ASCII')" << endl;
  fileSIMClusterM << "save('" << fileLSIMclusterStatistic << "','netsim','-ASCII','dpsim','-ASCII','expref','-ASCII')" << endl;
  fileSIMClusterM << "quit" << endl;
  fileSIMClusterM.close();

  sprintf(fileLSIMclusterName, "LSIMclusterTrial%03d.csh\0", posPart);

  //Execution of Script
  char *command = new char[1000];
  sprintf(command, "csh %s\0", fileLSIMclusterName); system(command);

  //Verification of the convergence process
  char *convergenceFile = new char[1000];
  sprintf(convergenceFile, "matlab%d.log\0", posPart);
  convergence = verifyAPAMatLab(convergenceFile);
  delete[] convergenceFile;

  if (!convergence) {
      
    sprintf(command, "rm LSIMclusterResult%03d\0",    posPart); system(command);
    sprintf(command, "rm LSIMclusterStatistic%03d\0", posPart); system(command);

    //Dinamic construction of matlab code
    sprintf(fileLSIMclusterName, "LSIMclusterTrial%03d.m\0", posPart);

    fileSIMClusterM.open(fileLSIMclusterName, fstream::out);

    fileSIMClusterM << "load " << similaritiesFolder << endl;
    if (USE_PREFEREN) {
      fileSIMClusterM << "load " << preferenceFolder << endl;
      fileSIMClusterM << "[idx,netsim,dpsim,expref]=apcluster(" << similaritiesFolder << "," << preferenceFolder << ",'maxits'," << MAXITS_2 << ",'convits'," << CONVITS_2 << ",'dampfact'," << DAMFACT_2 << ");" << endl;
    }
    else {
      fileSIMClusterM << "[idx,netsim,dpsim,expref]=apcluster(" << similaritiesFolder << "," << median << ",'maxits'," << MAXITS_2 << ",'convits'," << CONVITS_2 << ",'dampfact'," << DAMFACT_2 << ");" << endl;
    }
    fileSIMClusterM << "save('" << fileLSIMclusterResult << "','idx','-ASCII')" << endl;
    fileSIMClusterM << "save('" << fileLSIMclusterStatistic << "','netsim','-ASCII','dpsim','-ASCII','expref','-ASCII')" << endl;
    fileSIMClusterM << "quit" << endl;
    fileSIMClusterM.close();

    sprintf(fileLSIMclusterName, "LSIMclusterTrial%03d.csh\0", posPart);

    //Execution of Script
    sprintf(command, "csh %s\0", fileLSIMclusterName); system(command);
    
  }

  loadIdxAPAMatLab(fileLSIMclusterResult, IdxArray);

  sprintf(command, "rm matlab%d.log\0",             posPart); system(command);
  sprintf(command, "rm LSIMclusterTrial%03d.csh\0", posPart); system(command);
  sprintf(command, "rm LSIMclusterTrial%03d.m\0",   posPart); system(command);
  sprintf(command, "rm LSIMclusterResult%03d\0",    posPart); system(command);
  sprintf(command, "rm LSIMclusterStatistic%03d\0", posPart); system(command);

  delete[] fileLSIMclusterName;
  delete[] fileLSIMclusterResult;
  delete[] fileLSIMclusterStatistic;
  delete[] matlabPATH;
  delete[] command;

  return 1;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int computeAPA(char *similaritiesFolder, char *preferenceFolder, double median, IdxCluster *IdxArray, int length) {

  int firstExec   = true, iterations  = true, convergence = true;


  char *folder = new char[strlen(similaritiesFolder) + 1];
  strcpy(folder, similaritiesFolder);
  folder[strlen(folder) - 4] = '\0';

  while (iterations) {
    char *command = new char[1000];
    if (USE_PREFEREN) {
      if (firstExec) sprintf(command, "./apcluster %s %s %s maxits=%d convits=%d dampfact=%f\0", similaritiesFolder, preferenceFolder, folder, MAXITS_1, CONVITS_1, DAMFACT_1);
      else           sprintf(command, "./apcluster %s %s %s maxits=%d convits=%d dampfact=%f\0", similaritiesFolder, preferenceFolder, folder, MAXITS_2, CONVITS_2, DAMFACT_2);
    }
    else {
      if (USE_MEDIAN) {
        if (firstExec) sprintf(command, "./apcluster %s %f %s maxits=%d convits=%d dampfact=%f\0", similaritiesFolder, MEDIAN, folder, MAXITS_1, CONVITS_1, DAMFACT_1);
        else           sprintf(command, "./apcluster %s %f %s maxits=%d convits=%d dampfact=%f\0", similaritiesFolder, MEDIAN, folder, MAXITS_2, CONVITS_2, DAMFACT_2);
      }
      else {
        if (firstExec) sprintf(command, "./apcluster %s %f %s maxits=%d convits=%d dampfact=%f\0", similaritiesFolder, median, folder, MAXITS_1, CONVITS_1, DAMFACT_1);
        else           sprintf(command, "./apcluster %s %f %s maxits=%d convits=%d dampfact=%f\0", similaritiesFolder, median, folder, MAXITS_2, CONVITS_2, DAMFACT_2);
      }
    }
    cout << command << endl;
    system(command);
    if (command) delete[] command;

    convergence = verifyAPA(folder);
    if ((!convergence) || (convergence && !firstExec)) break;
    firstExec = false;
  }

  if (!convergence)
    loadIdxAPA(folder, IdxArray);
  else
    for (int i = 0; i < length; i++)
      IdxArray[i].IdxValue(i + 1);

  char *command = new char[1000];
  sprintf(command, "rm %s/*.*;rmdir %s\0", folder, folder);
  system(command);
  if (command) delete[] command;
  if (folder) delete[] folder;

  return 1;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int clusteringProcess(smileStore *molDB, char *filename, int totalNumberSmiles, int sizeRank, int rank, char *resultFile) {

  int initialPos, finalPos;                                                           //Initial and Final Position of Smiles in the original File
  int numberPart = totalNumberSmiles/CLUST_SIZE;                                      //Number of Data Partition
  int iter = 0, posPart = 0;                                                          //Number of Iteration to complete the operations and Partition Position

  if (!numberPart && rank) return 0;

  numberPart +=  (totalNumberSmiles%CLUST_SIZE)?1:0;

  molDB->Filename(filename);

  posPart = ((iter++)*sizeRank) + rank;                                               //Partition Position

  if (posPart > numberPart - 1) return 0;

  //echo '	The complete '$DB.smi' database containing '$TOTcompN' entries is divided'
  //echo '	into '$RGCgroupsN' random groups of '$ClustGroupsCapacity' compounds'

  //Similiraties File
  fstream fileResult;
  sprintf(resultFile, "tmp_Result%0.5d.txt\0", rank);
  fileResult.open(resultFile, fstream::in | fstream::out | fstream::app);

  while (posPart < numberPart) {

    initialPos = posPart*CLUST_SIZE;
    finalPos   = ((posPart + 1) < numberPart)?((posPart + 1)*CLUST_SIZE - 1):(totalNumberSmiles - 1);

    molDB->ExtractSmiles(initialPos, finalPos);                                       //Extract only the smiles to be processing in this iterations

    //Similiraties File
    fstream fileout;
    char    *similaritiesFile = new char[30];
    char    *preferenceFile   = new char[30];
    sprintf(similaritiesFile, "Similarities%0.5d.txt\0", posPart);
    sprintf(preferenceFile, "Preference%0.5d.txt\0", posPart);
    fileout.open(similaritiesFile, ios::out);

    char *buffer     = new char[LENGHTBLINE*molDB->NumberSmiles() + 1];
    char *bufferLine = new char[LENGHTBLINE + 1];
    buffer[0] = '\0';

    //Preference File
    if (USE_PREFEREN) {
      fstream fileoutP;
      fileoutP.open(preferenceFile, ios::out);

      for (int i = 0; i < molDB->NumberSmiles(); i++) {
        if (!i) sprintf(bufferLine, "%.6f", MEDIAN);
        else    sprintf(bufferLine, "\n%.6f", MEDIAN);
        sprintf(buffer, "%s%s\0", buffer, bufferLine);
      }

      fileoutP << buffer;
      fileoutP.close();
    }
    
    if (buffer) delete[] buffer;
    buffer    = new char[LENGHTBLINE*BUFFCICLE + 1];
    buffer[0] = '\0';

    int     counterCicles = 0;
    double  medianValue   = 0;

    //Compute the LingoSim for the group of smiles in DB
    for (int i = 0; i < molDB->NumberSmiles(); i++) 
      for (int j = 0; j < i; j++) {

        //double tmp_value  = ComputeLingoSimilarity(molDB->OneSmileStruct(i), molDB->OneSmileStruct(j));
        double tmp_value = obtainLingoSim(molDB->OneSmileStruct(i)->Code(), molDB->OneSmileStruct(j)->Code());

        if (USE_DISTANCE) tmp_value = tmp_value - 1;
        
        medianValue += tmp_value;

        sprintf(bufferLine, "%04d  %04d  %.6f\n", i + 1, j + 1, tmp_value);
        sprintf(buffer, "%s%s\0", buffer, bufferLine);

        sprintf(bufferLine, "%04d  %04d  %.6f\n", j + 1, i + 1, tmp_value);
        sprintf(buffer, "%s%s\0", buffer, bufferLine);

        counterCicles += 2;

        if (counterCicles >= BUFFCICLE) {
          counterCicles = 0;
          fileout << buffer;
          buffer[0] = '\0';
        }

      }

    medianValue *= 2;
    medianValue += (USE_DISTANCE)?0:molDB->NumberSmiles();
    medianValue /= (molDB->NumberSmiles()*molDB->NumberSmiles());

    if (counterCicles < BUFFCICLE) {
      counterCicles = 0;
      fileout << buffer;
      buffer[0] = '\0';
    }
      
    if (buffer) delete[] buffer;
    
    fileout.close();
    
    //Compute the clusters with Affinity Propagation Algorimths
    IdxCluster *IdxArray = new IdxCluster[molDB->NumberSmiles()];
    for (int i = 0; i < molDB->NumberSmiles(); i++) IdxArray[i].IdxValue(-1);

    computeAPA(similaritiesFile, preferenceFile, medianValue, IdxArray, molDB->NumberSmiles());

    //Load the exemplars and store the trace

    for (int i = 0; i < molDB->NumberSmiles(); i++) {
      fileResult << left;
      fileResult.width(12);
      fileResult << molDB->OneSmileStruct(i)->Id() << "  " << molDB->OneSmileStruct(IdxArray[i].IdxValue() - 1)->Id() << endl;
    }

    char *cmdDelFile = new char[strlen(similaritiesFile) + 10];
    sprintf(cmdDelFile, "rm %s\0", similaritiesFile);
    system(cmdDelFile);
    if (USE_PREFEREN) {
      sprintf(cmdDelFile, "rm %s\0", preferenceFile);
      system(cmdDelFile);
    }
    if (similaritiesFile) delete[] similaritiesFile;
    if (preferenceFile) delete[] preferenceFile;

    if (IdxArray) delete[] IdxArray;
    
    posPart = ((iter++)*sizeRank) + rank;                                                  //Partition Position
  

  }
  fileResult.close();

  return 1;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int clusteringProcessMatLab(smileStore *molDB, char *filename, int totalNumberSmiles, int sizeRank, int rank, char *resultFile) {

  int initialPos, finalPos;                                                           //Initial and Final Position of Smiles in the original File
  int numberPart = totalNumberSmiles/CLUST_SIZE;                                      //Number of Data Partition
  int iter = 0, posPart = 0;                                                          //Number of Iteration to complete the operations and Partition Position

  if (!numberPart && rank) return 0;

  numberPart +=  (totalNumberSmiles%CLUST_SIZE)?1:0;

  molDB->Filename(filename);

  posPart = ((iter++)*sizeRank) + rank;                                               //Partition Position

  if (posPart > numberPart - 1) return 0;

  //Similiraties File
  fstream fileResult;
  sprintf(resultFile, "tmp_Result%0.5d.txt\0", rank);
  fileResult.open(resultFile, fstream::in | fstream::out | fstream::app);

  while (posPart < numberPart) {

    initialPos = posPart*CLUST_SIZE;
    finalPos   = ((posPart + 1) < numberPart)?((posPart + 1)*CLUST_SIZE - 1):(totalNumberSmiles - 1);

    molDB->ExtractSmiles(initialPos, finalPos);                                       //Extract only the smiles to be processing in this iterations

    //Similiraties File
    fstream fileout;
    char    *similaritiesFile = new char[30];
    char    *preferenceFile   = new char[30];
    sprintf(similaritiesFile, "Similarities%0.5d\0", posPart);
    sprintf(preferenceFile, "Preference%0.5d\0", posPart);
    fileout.open(similaritiesFile, ios::out);

    char *buffer     = new char[LENGHTBLINE*molDB->NumberSmiles() + 1];
    char *bufferLine = new char[LENGHTBLINE + 1];
    buffer[0] = '\0';

    //Preference File
    if (USE_PREFEREN) {
      fstream fileoutP;
      fileoutP.open(preferenceFile, ios::out);

      for (int i = 0; i < molDB->NumberSmiles(); i++) {
        if (!i) sprintf(bufferLine, "%.6f", MEDIAN);
        else    sprintf(bufferLine, "\n%.6f", MEDIAN);
        sprintf(buffer, "%s%s\0", buffer, bufferLine);
      }

      fileoutP << buffer;
      fileoutP.close();
    }

    if (buffer) delete[] buffer;
    buffer    = new char[LENGHTBLINE*BUFFCICLE + 1];
    buffer[0] = '\0';

    int     counterCicles = 0;
    double  medianValue   = 0;

    //Compute the LingoSim for the group of smiles in DB
    for (int i = 0; i < molDB->NumberSmiles(); i++)
      for (int j = 0; j < i; j++) {

        //double tmp_value  = ComputeLingoSimilarity(molDB->OneSmileStruct(i), molDB->OneSmileStruct(j));
        double tmp_value = obtainLingoSim(molDB->OneSmileStruct(i)->Code(), molDB->OneSmileStruct(j)->Code());

        if (USE_DISTANCE) tmp_value = tmp_value - 1;

        medianValue += tmp_value;

        sprintf(bufferLine, "%04d  %04d  %.6f\n", i + 1, j + 1, tmp_value);
        sprintf(buffer, "%s%s\0", buffer, bufferLine);

        sprintf(bufferLine, "%04d  %04d  %.6f\n", j + 1, i + 1, tmp_value);
        sprintf(buffer, "%s%s\0", buffer, bufferLine);

        counterCicles += 2;

        if (counterCicles >= BUFFCICLE) {
          counterCicles = 0;
          fileout << buffer;
          buffer[0] = '\0';
        }

      }

    medianValue *= 2;
    medianValue += (USE_DISTANCE)?0:molDB->NumberSmiles();
    medianValue /= (molDB->NumberSmiles()*molDB->NumberSmiles());

    if (counterCicles < BUFFCICLE) {
      counterCicles = 0;
      fileout << buffer;
      buffer[0] = '\0';
    }

    if (buffer) delete[] buffer;
    fileout.close();

    //Compute the clusters with Affinity Propagation Algorimths
    IdxCluster *IdxArray = new IdxCluster[molDB->NumberSmiles()];
    for (int i = 0; i < molDB->NumberSmiles(); i++) IdxArray[i].IdxValue(-1);

    computeAPAMatLab(similaritiesFile, preferenceFile, medianValue, IdxArray, molDB->NumberSmiles(), posPart);

    //Load the exemplars and store the trace
    for (int i = 0; i < molDB->NumberSmiles(); i++) {
      if (IdxArray[i].IdxValue() == -1) continue;
      fileResult << left;
      fileResult.width(12);
      fileResult << molDB->OneSmileStruct(i)->Id() << "  " << molDB->OneSmileStruct(IdxArray[i].IdxValue() - 1)->Id() << endl;
    }

    char *cmdDelFile = new char[strlen(similaritiesFile) + 10];
    sprintf(cmdDelFile, "rm %s\0", similaritiesFile);
    system(cmdDelFile);
    if (USE_PREFEREN) {
      sprintf(cmdDelFile, "rm %s\0", preferenceFile);
      system(cmdDelFile);
    }
    
    if (similaritiesFile) delete[] similaritiesFile;
    if (preferenceFile) delete[] preferenceFile;

    if (IdxArray) delete[] IdxArray;

    posPart = ((iter++)*sizeRank) + rank;                                                  //Partition Position


  }
  fileResult.close();

  return 1;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int splitLineIdCode(char *line, int *idCodes, int size) {

  char *pch, *lineTmp = new char[strlen(line) + 1];
  int   count = 0;
  
  strcpy(lineTmp, line);

  pch              = strtok(lineTmp, " ");
  idCodes[count++] = atoi(pch);

  while (pch != NULL) {
    pch              = strtok(NULL, " ");
    if (count < size) idCodes[count++] = atoi(pch);
  }

  delete[] lineTmp;

  return 1;
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int splitPosLineIdCode(char *line, int pos) {

  char *pch, *lineTmp = new char[strlen(line) + 1];
  int   count = 1;

  strcpy(lineTmp, line);

  pch = strtok(lineTmp, " ");
  if (pos == 1) return atoi(pch);

  while (pch != NULL) {
    pch = strtok(NULL, " ");
    if (++count == pos) return atoi(pch);
  }

  delete[] lineTmp;

  return -1;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int joinClusterResult(char *inputNameFile, int nAncestor, int totalSmiles, int newTotalSmiles) {

  if (!inputNameFile || nAncestor < 0) return 0;

  if (!nAncestor) {
    //Create backup
    char *command = new char[100];
    sprintf(command, "cp %s Cluster_Molecules.txt\0", inputNameFile); system(command);
    if (command) delete[] command;
  }
  else if (nAncestor > 0) {

    fstream clusterFile, inputFile;
    clusterFile.open("Cluster_Molecules.txt", fstream::in | fstream::out | fstream::app);
    inputFile.open(inputNameFile, fstream::in);

    int **lingosId = new int*[totalSmiles];
    for (int i = 0; i < totalSmiles; i++) lingosId[i] = new int[nAncestor + 2];

    int **inputId = new int*[newTotalSmiles];
    for (int i = 0; i < newTotalSmiles; i++) inputId[i] = new int[2];

    //Read Cluster_Molecules
    clusterFile.seekg (0, fstream::beg);
    int count = 0;
    char *line = new char[1000];
    while (!clusterFile.eof()) {
      clusterFile.getline(line, 1000);
      if (!strlen(line)) continue;
      splitLineIdCode(line, lingosId[count++], nAncestor + 1);
    }

    //Read Input_File
    inputFile.seekg (0, fstream::beg);
    count = 0;
    while (!inputFile.eof()) {
      inputFile.getline(line, 1000);
      if (!strlen(line)) continue;
      splitLineIdCode(line, inputId[count++], 2);
    }
    if (line) delete[] line;

    //Complete Cluster_Molecule
    for (int i = 0; i < totalSmiles; i++) {
      bool condition = false;
      for (int j = 0; j < newTotalSmiles; j++)
        if (lingosId[i][nAncestor] == inputId[j][0]) {
          lingosId[i][nAncestor + 1] = inputId[j][1];
          condition = true;
          break;
        }
    }

    clusterFile.close();
    clusterFile.open("Cluster_Molecules.txt", fstream::out);

    for (int i = 0; i < totalSmiles; i++) {
      for (int j = 0; j < (nAncestor + 2); j++) {
        clusterFile << left;
        clusterFile.width(12);
        clusterFile << lingosId[i][j] << " ";
      }
      clusterFile << endl;
    }

    for (int i = 0; i < totalSmiles; i++) if (lingosId[i]) delete[] lingosId[i];
    if (lingosId) delete[] lingosId;

    for (int i = 0; i < newTotalSmiles; i++) if (inputId[i]) delete[] inputId[i];
    if (inputId) delete[] inputId;
    
    inputFile.close();
    clusterFile.close();

  }
    
  return 1;
    
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int prepareNewDB(char *dataFile, char *outFile, int totalNumberSmiles, smileStore *molDB) {

  fstream fileResult;
  fileResult.open(dataFile, fstream::in);
  fileResult.seekg (0, fstream::beg);

  int *codeId = new int[totalNumberSmiles];

  molDB->Filename(outFile);
  molDB->ExtractSmiles();

  int   count = 0;
  char *line = new char[1000];
  while (!fileResult.eof()) {
    fileResult.getline(line, 1000);
    if (!strlen(line)) continue;
    codeId[count++] = splitPosLineIdCode(line, 2);
  }
  if (line) delete[] line;

  smileStruct::MergeSort(codeId, totalNumberSmiles);
  totalNumberSmiles = smileStruct::CompactList(codeId, totalNumberSmiles);
  fileResult.close();

  fstream newFile;
  sprintf(outFile, "DB_Molecule_TMP.tmp\0");
  newFile.open(outFile, fstream::out);

  for (int i = 0; i < totalNumberSmiles; i++)
    for (int j = 0; j < molDB->NumberSmiles(); j++)
      if (atoi(molDB->OneSmileStruct(j)->Id()) == codeId[i]) {
          newFile << molDB->OneSmileStruct(j)->Code() << " " << molDB->OneSmileStruct(j)->Id() << endl;
          break;
      }

  newFile.close();

  if (codeId) delete[] codeId;

  return 1;
    
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int loadMoleculeCluster(char *filename, long **&listClusterFile, int numberMolecules) {

  if (!filename) return 0;
  
  fstream fileInput;
  fileInput.open(filename, fstream::in);

  int   count = 0;
  
  while (!fileInput.eof()) {
    char *line  = new char[1001];
    fileInput.getline(line, 1000);
    if (strlen(line)) {
      listClusterFile[count][0]   = splitPosLineIdCode(line, 1);
      listClusterFile[count++][1] = splitPosLineIdCode(line, 2);
    }
    if (count > numberMolecules) break;
    delete[] line;
  }
  
  
  fileInput.close();

  return 1;
    
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int findCommonData(long *list1, long *list2, int nlist1, int nlist2) {

  int numberCommonData = 0;
  for (int i = 0; i < nlist1; i++) for (int j = 0; j < nlist2; j++) if (list1[i] == list2[j]) { numberCommonData++; break; }
  return numberCommonData;
    
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double computeConvergence(long **&listClusterFileOld, long **&listClusterFileNew, int numberMolecules) {

  long   *listNew = new long[numberMolecules];
  long   *listOld = new long[numberMolecules];
  
  double convergenceValue = 0;
  double factor, fnumberCommonCluster, fnumberOldCluster, fnumberNewCluster;

  for (int i = 0; i < numberMolecules; i++) {
      
    //Initializate the subList of Molecule
    for (int j = 0; j < numberMolecules; j++) listOld[j] = listNew[j] = -1;
    
    //Get the father fot the New Cluster of Molecules
    int numberNewCluster = 0, numberOldCluster = 0, numberCommonCluster = 0;
    
    for (int j = 0; j < numberMolecules; j++) 
      if (listClusterFileNew[i][1] == listClusterFileNew[j][1]) 
        listNew[numberNewCluster++] = listClusterFileNew[j][0];
    
    long fatherMolOld = -1;
    
    for (int j = 0; j < numberMolecules; j++)
      if (listClusterFileNew[i][0] == listClusterFileOld[j][0]) {
        fatherMolOld = listClusterFileOld[j][1];
        break;
      }
    
    for (int j = 0; j < numberMolecules; j++)
      if (fatherMolOld == listClusterFileOld[j][1])
        listOld[numberOldCluster++] = listClusterFileOld[j][0];

    if (numberOldCluster < numberNewCluster) numberCommonCluster = findCommonData(listOld, listNew, numberOldCluster, numberNewCluster);
    else                                     numberCommonCluster = findCommonData(listNew, listOld, numberNewCluster, numberOldCluster);

    fnumberCommonCluster = numberCommonCluster; fnumberOldCluster = numberOldCluster; fnumberNewCluster = numberNewCluster;
    factor = (fnumberCommonCluster - 1)/(fnumberOldCluster + fnumberNewCluster - fnumberCommonCluster);
    convergenceValue += factor;

  }

  delete[] listNew;
  delete[] listOld;

  return convergenceValue;
    
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool analizeConvergence(int iteration, int nSmiles, int rank, double &convergenceValue) {

  if (!iteration) {
    if (!rank) {
      char *command = new char[100];
      sprintf(command, "cp Cluster_Molecules.txt Cluster_Molecules%03d.txt\0", iteration); system(command);
      sprintf(command, "rm Cluster_Molecules.txt\0"); system(command);
      if (command) delete[] command;
    }
    return false;
  }

  if (!rank) {

    static int count = 0;
    long **listClusterFileOld = new long*[nSmiles];
    long **listClusterFileNew = new long*[nSmiles];

    for (int i = 0; i < nSmiles; i++) {
      listClusterFileOld[i] = new long[2];
      listClusterFileNew[i] = new long[2];
    }

    char *command = new char[100];
    sprintf(command, "cp Cluster_Molecules.txt Cluster_Molecules%03d.txt\0", iteration); system(command);
    sprintf(command, "rm Cluster_Molecules.txt\0"); system(command);

    sprintf(command, "Cluster_Molecules%03d.txt\0", iteration - 1);
    loadMoleculeCluster(command, listClusterFileOld, nSmiles);

    sprintf(command, "Cluster_Molecules%03d.txt\0", iteration);
    loadMoleculeCluster(command, listClusterFileNew, nSmiles);

    if (command) delete[] command;

    double valueConvergence = computeConvergence(listClusterFileOld, listClusterFileNew, nSmiles);
    double ownConvergence   = computeConvergence(listClusterFileNew, listClusterFileNew, nSmiles);

    for (int i = 0; i < nSmiles; i++) {
      if (listClusterFileOld[i]) delete[] listClusterFileOld[i];
      if (listClusterFileNew[i]) delete[] listClusterFileNew[i];
    }
    if (listClusterFileOld) delete[] listClusterFileOld;
    if (listClusterFileNew) delete[] listClusterFileNew;

    double percentConvergence = PERCENTCONV;
    convergenceValue = valueConvergence/ownConvergence;
    if (convergenceValue >= percentConvergence) return true;
    return false;

  }
  return false;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int main (int argc, char *argv[]) {

  time_t initialTime, finalTime;

  initialTime = time(NULL);

  int           err, rank, sizeRank;            //Variable to use in MPI functions
  MPI_Request	reqs[ARRAYSIZE];	        //Variable to use in MPI functions
  MPI_Status	Stat;                           //Variable to use in MPI functions

  smileStore	*molDB = new smileStore;	//Variable to store the smiles and lingo datas for all molecules in DB

  char 	  filename[100], initialFileName[100];
  char    nameLogFile[] = {"parallel_exec_LINGOSim.log"};
  int     lengthFileName, totalNumberSmiles, nSmiles;
  int     nconvergence = 0;
  bool    convergence  = false, condLogFile = false;

  fstream logFile;

  err = MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &sizeRank);

  if (!rank) {
    logFile.open(nameLogFile, fstream::out);
    condLogFile = logFile.is_open();
    if (condLogFile) logFile << "Clustering Process Starting....." << endl << "Number of Parallel Execution Process (" << sizeRank << ")" << endl;
  }

  //Iteration until the clustering process converg
  while (!convergence) {

    if (nconvergence && !rank) strcpy(filename, initialFileName);
    MPI_Barrier(MPI_COMM_WORLD);

    if (!rank && condLogFile) logFile << "Iteration: " << nconvergence << endl;

    //Analysis of ancestor for each convegence process
    for (int ancestor = 0; ancestor < NANCESTOR; ancestor++) {

      //Master Process
      if (!rank) {
          
        if (condLogFile) logFile << "Ancestral Level: " << ancestor << endl;
        
        //Input data for the first time of execution
        if (!ancestor && !nconvergence) {
          if      (argc < 2)  { printf("Enter the input file: \t"); gets(filename); }
          else if (argc == 2) { strcpy(filename, argv[1]); }
          printf("The DB file is: %s\n", filename);
          if (condLogFile) logFile << "The DB file is: " << filename << endl;
          strcpy(initialFileName, filename);
        }

        //Send the filename variable to all workers process
        lengthFileName = strlen(filename);
        
        for (int i = 1; i < sizeRank; i++) {
          MPI_Send(&lengthFileName,              1,  MPI_INT, i, (10 + i), MPI_COMM_WORLD);
          MPI_Send(filename       , lengthFileName, MPI_CHAR, i, (20 + i), MPI_COMM_WORLD);
        }
        
      }
      //Worker Process
      else {

        //Receive the filename variable
        MPI_Recv(&lengthFileName,	       1,  MPI_INT, 0, (10 + rank), MPI_COMM_WORLD, &Stat);
        MPI_Recv(filename       , lengthFileName, MPI_CHAR, 0, (20 + rank), MPI_COMM_WORLD, &Stat);
        filename[lengthFileName] = '\0';
        molDB->Filename(filename);

      }

      molDB->Filename(filename);

      //Prepare the input File (Only in the Master Process)
      if (!rank) {

        if (condLogFile) logFile << "Preparation of Input File....." << endl;

        if (!nconvergence) {                                                             //First execution of algorithm for all ancestral level (Ramdon Sequences Order)
          if (condLogFile) logFile << "Order the sequences in Ramdon Order....." << endl;
          molDB->RandomFile(ITERATIONS);
        }
        else if (nconvergence && !ancestor) {                                           //Execution of algorithm for the first ancestral level after the first execution (Sequences Ordering by level of ancestral)
          if (condLogFile) logFile << "Order the sequences by the Previos Clustering....." << endl;
          char *nameFileClusterMol = new char[100];
          sprintf(nameFileClusterMol, "Cluster_Molecules%03d.txt\0", nconvergence - 1);
          molDB->OrderFile(filename, nameFileClusterMol, NANCESTOR);
          delete[] nameFileClusterMol;
        }
        else {                                                                          //Execution of algorithm fot the level different of first ancestral level after the first execution (Sequences Ordering keeping the order of previous levels)
          if (condLogFile) logFile << "Order the sequences keeping the order in the Previos Clustering....." << endl;
          char *nameFileClusterMol = new char[100];
          sprintf(nameFileClusterMol, "Cluster_Molecules.txt\0");
          molDB->ReduceFile(filename, nameFileClusterMol);
          delete[] nameFileClusterMol;
        }

        totalNumberSmiles = molDB->NumberSmiles();

        for (int i = 1; i < sizeRank; i++)
          MPI_Send(&totalNumberSmiles, 1,  MPI_INT, i, (10 + i), MPI_COMM_WORLD);
        
        if (!ancestor) nSmiles = totalNumberSmiles;
        
      }
      else
        MPI_Recv(&totalNumberSmiles, 1,  MPI_INT, 0, (10 + rank), MPI_COMM_WORLD, &Stat);

      //Clustering process
      if (!rank && condLogFile) logFile << "Clustering Process Starting....." << endl;
      char *resultFile = new char[30];
      int   stateWork;
      if (USE_MATLAB) stateWork  = clusteringProcessMatLab(molDB, filename, totalNumberSmiles, sizeRank, rank, resultFile);
      else            stateWork  = clusteringProcess(molDB, filename, totalNumberSmiles, sizeRank, rank, resultFile);
      if (!rank && condLogFile) logFile << "Clustering Process Finishing....." << endl;

      if (!rank && condLogFile) logFile << "Synchronization Phase....." << endl;
      MPI_Barrier(MPI_COMM_WORLD);

      if (!rank && condLogFile) logFile << "Communication Phase....." << endl;
      //Integrate the obtained result for each node (if the node work with any data)
      if (stateWork) {

        if (!rank) {

          fstream fileResult;
          sprintf(resultFile, "tmp_Result%0.5d.txt\0", rank);
          fileResult.open(resultFile, fstream::in | fstream::out | fstream::app);
          fileResult.seekg (0, fstream::end);

          long  length = 0;
          char *buffer = NULL;

          for (int i = 1; i < sizeRank; i++) {
            MPI_Recv(&length,      1,  MPI_INT, i, (30 + i), MPI_COMM_WORLD, &Stat);
            buffer = new char[length + 1];
            MPI_Recv(buffer , length, MPI_CHAR, i, (40 + i), MPI_COMM_WORLD, &Stat);
            buffer[length] = '\0';
            fileResult << buffer;
            if (buffer) delete[] buffer;
          }

          fileResult.close();

          char *command = new char[100];

          if (condLogFile) logFile << "Joining the Cluster Result....." << endl;
          joinClusterResult(resultFile, ancestor, nSmiles, totalNumberSmiles);
          if (condLogFile) logFile << "Preparing the new Data Set....." << endl;
          prepareNewDB(resultFile, filename, totalNumberSmiles, molDB);

          sprintf(command, "rm %s\0", resultFile); system(command);
          if (command) delete[] command;

        }
        else {

          fstream fileResult;
          sprintf(resultFile, "tmp_Result%0.5d.txt\0", rank);
          fileResult.open(resultFile, fstream::in);

          fileResult.seekg (0, fstream::end);
          long length = fileResult.tellg();
          fileResult.seekg (0, fstream::beg);
          char *buffer = new char[length + 1];

          //Read data as a block:
          fileResult.read(buffer, length);
          fileResult.close();

          MPI_Send(&length,      1,  MPI_INT, 0, (30 + rank), MPI_COMM_WORLD);
          MPI_Send(buffer , length, MPI_CHAR, 0, (40 + rank), MPI_COMM_WORLD);

          if (buffer) delete[] buffer;

          char *command = new char[40];
          sprintf(command, "rm %s\0", resultFile); system(command);
          if (command) delete[] command;

        }

      }
      else {

        long  length = 0;
        char *buffer = new char[1];

        MPI_Send(&length,      1,  MPI_INT, 0, (30 + rank), MPI_COMM_WORLD);
        MPI_Send(buffer , length, MPI_CHAR, 0, (40 + rank), MPI_COMM_WORLD);

        if (buffer) delete[] buffer;

      }

      if (resultFile) delete[] resultFile;

      MPI_Barrier(MPI_COMM_WORLD);

    }
      
    if (!rank) {

      char *command = new char[100];
      sprintf(command, "rm %s\0", filename);
      system(command);
      if (command) delete[] command;
      
    }

    //Analysis the convergence with the previous iteration
    if (!rank && condLogFile) logFile << "Analysis the Convergence with the Previous Iteration....." << endl;
    double convergenceValue = 0;
    convergence = analizeConvergence(nconvergence, nSmiles, rank, convergenceValue);
    if (!rank && condLogFile) logFile << "Convergence = " << convergenceValue << "(Iteration: " << nconvergence << ")" << endl;
    
    nconvergence++;

  }
  
  MPI_Finalize();

  if (condLogFile) logFile << "Parallel Process Finishing....." << endl;
    
  if (molDB) delete molDB;

  finalTime = time(NULL);
  if (condLogFile) {
    logFile << "Time in Seconds: (" << finalTime - initialTime << ")" << endl;
    logFile.close();
  }
  
  return 1;
 	
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
